Project supported by the National Key Research and Development Program of China (Grant No. 2017YFB1104500), the Natural Science Foundation of Beijing (Grant No. 7182091), and the National Natural Science Foundation of China (Grant No. 21627813).
Abstract
This paper puts forward for the first time a combined transmission matrix (TM) method to measure the monochromatic TM of scattering media without a reference beam. This method can be named a sequential semi-definite programming method which combines the sequential algorithm and the semi-definite programming method. Firstly, each part of the TM is calculated respectively in proper sequence. Then every part of TM is combined to form a complete TM in accordance with a certain rule. The phase modulation of the incident light is achieved by using a high speed digital mirror device with the superpixel method. We have experimentally demonstrated that the incident light field is focused at the target through scattering media using the measured TM to optimize the wavefront of the incident light. Compared with the semi-definite programming method, our method takes less computational time and occupies less memory space. The sequential semi-definite programming method shows potential applications in imaging through biological tissues.
Wave propagation in complex media is a fundamental problem in acoustic, electromagnetism, and optics.[1] Transmission matrix (TM) is the basic model to describe wave propagation within a medium. In optics, controlling light propagation in scattering media has a wide range of application, such as biological imaging, information security, and ghost imaging.[2–4] In 2007, Vellekoop et al. put forward a method of wavefront-shaping for the first time. They focused the incident beam through scattering media with iterative optimization algorithm.[5] There are many algorithms to optimize the wavefront, such as genetic algorithm, sequential algorithm, and particle swarm optimization algorithm.[6–8] All of these algorithms require feedback to optimize the wavefront. In 2010, Popoff et al. proposed another method for wavefront-shaping.[9] They reconstructed an original image from the speckle field by measuring the monochromatic TM of multiple scattering media. In this TM measurement method, feedback regulation was not a necessity. However, the traditional TM measurement methods, such as four phase method[10] and Lee holographic method,[11] needed to measure the output complex field by interfering with the reference light field, which made the experimental setup more complex. In recent years, many phase retrieval algorithms were applied to measure TM of multiple scattering media, such as prGAMP, prSAMP, prVBEM, prVAMP, GS, WF, PhaseLift, and PhaseMax.[12] But the author used a high-throughput computing (HTC) cluster to compute TM of multiple scattering that personal computer could not complete. In 2017, N’Gom et al. measured the TM with the semi-definite programming (SDP) method,[13] which simplified the experimental setup because there was no reference beam and personal computer could complete the work. However, with the increase of TM elements, the method required more calculation time and consumed more memory space. For example, with the element number of 36, the method took 10.4 minutes and consumed 2-GB memory space. However, with the element number of 49, the method took 106 minutes and consumed 40 GB memory space. In consequence, as the number of transmission matrix elements increased, the calculation time and the memory space increased exponentially. Therefore, the computational complexity of this method was a great challenge for personal computers when the number of elements was large.
In this paper, we propose for the first time a combined transmission matrix method to measure the monochromatic TM of scattering media without a reference beam, i.e., the sequential semi-definite programming (SSDP) method. In this SSDP method, we first divide the TM into multiple parts and calculate each part in proper sequence. Then we assemble all parts together and obtain a complete TM after a corrective optimization. Therefore, the SSDP method can greatly reduce the calculation time and memory. This method also provides a new experimental idea to measure large TM of scattering medium because it can reduce the calculation time and memory. In the current experiments, we use digital micromirror device (DMD) to modulate the phase of the incident wavefront with the superpixel method.[14,15] The SSDP method is used to compute the TM of the scattering medium by recording only the intensity information of the output light field. Then based on this TM, we compute the optimal wavefront for focusing the incident light through the scattering medium. Our method has greatly improved the speed of light focusing compared with the SDP method.
2. Method
2.1. Superpixel method
We use DMD as spatial light modulator, which is an excellent electronic equipment for controlling light fields. It owns millions of switchable micromirror units and has an extremely fast refresh rate up to 23 kHz. However, the DMD can only achieve spatial light binary amplitude modulation in general and we need phase modulation to get random phase mask in the SSDP method. Thus, we use the superpixel method to achieve spatial light phase modulation on the DMD.[14,15] The superpixel method, which is simple and highly robust, takes full advantage of the characteristics of DMD. In this method, the adjacent m × m micromirrors are grouped into single surperpixel and each superpixel defines a complex field in the target plane. After a plane light beam illuminates DMD surface with an angle, the reflected light is divided into many diffraction orders. In the +1 or −1 order, the phase is modulated by superpixel-based DMD. A 4 × 4 superpixel is shown in Fig. 1(a). If we turn on the three micromirrors indicated by green squares in Fig. 1(a), the corresponding superpixel in the target plane is the sum of three active pixels in Fig. 1(b). Using this method, we can achieve amplitude and phase modulation. In Fig. 1(c), we show a DMD pattern consisting of 9 × 9 segments with random phases. One segment is composed of many identical superpixels. In this paper, we just achieve phase-only modulation by the superpixel method.
Fig. 1. (color online) (a) Demonstration of a superpixel composed of 4 × 4 DMD mirrors which divide 2π into 16 parts. (b) The contribution of each individual mirror is located on a circle in the complex plane. (c) A DMD pattern consisting of 9 × 9 segments with random phase and one segment is composed of many identical superpixels.
2.2. SSDP method
The process of TM measurement by the SSDP method is shown in Fig. 2. A highly scattering medium is sandwiched between input field A and output field B. To clearly explain this method, we assume that both input field and output field consist of N = 16 × 16 segments. One input segment is composed of many identical superpixels. Light in the input field is characterized by N complex-valued electric field vectors an. Depending on the size of one superpixel, we can control the phase of an under different precisions. In the same way, light in the output field is characterized by N complex-valued electric vectors bn. The intensity of each bn can be recorded by CCD camera. Therefore, we need to retrieve the phase of bn.
Fig. 2. (color online) The process of TM measurement by the SSDP method. A is the input field, B is the output field, K is the scattering medium, an is one segment in the input field, and bn is one segment in the output field.
The input field and the output field are related by the N × N TM of the scattering medium
where A and B are N × 1 vectors of input field and output field. In order to get a focus on the output field, the TM of the scattering medium must be known. We assume that the light is focused at point bn, and the n-th row of the corresponding TM is Tn. Then equation (1) becomes
To determine Tn, we illuminate the scattering medium with M randomly selected fields A(m) (m = 1,…, M) of uniform magnitude. To get accurate Tn, the number of randomly selected fields must be M > N log (N), where N is the number of segments of the input field.
As the number of transmission matrix elements increases, the calculation time and the memory space increase exponentially. Therefore, computing costs much time and memory space when the number of segment of input field is large. We use the following method to reduce the calculation time and memory space. First, we divide Tn into several parts. In order to minimize the calculation time and memory space, the number of parts is always the square root of the number of Tn elements. If the square root is not an integer, we choose an integer that is just greater than the value of the square root. Second, we calculate each part of Tn separately. Third, every part of Tn combines to form a complete TM. As shown in Fig. 2, the Tn is divided into sixteen parts and each part has sixteen segments, such as Tn1 = [t1,t2,…,t16], Tn2 = [t17,t18,…,t32],…, Tn16 = [t241,t242,…,t256].
2.2.1. Single calculation of a part of TM
As Step 1 shown in Fig. 2, the magnitude of the field at the target is recorded for all excitations
and represents the phase information of . Then, m randomly selected fields (, i = 1, 2,…, m) are selected as input. The number of randomly selected fields is m > 16log (16). Next, we can get
where , , and . And the U1 and Tn1 can be obtained by solving the following optimization problem:
The above problem can be solved by using the SDP method[13] and the PhaseCut algorithm.[16] Equation (5) becomes:
where † is the Moor–Penrose pseudoinverse symbol, * is conjugate transpose symbol, and . Equation (6) is a convex optimization problem and it can be solved with numerical solvers by the cvx package. Then we can compute the phase information of target
where v1(V1) is the eigenvector of V1 associated with its largest eigenvalue. From Eq. (4), we can compute the TM
We next compute the optimization incident field vector A1 for the first part. The optimization input pattern can be written as[13]
2.2.2. Calculation of other parts of TM in sequence
For Step 2–Step 16 shown in Fig. 2, we follow Step 1 to compute other parts of transmission matrix Tn2,Tn3,…,Tn16, and compute the optimal incident field .
2.2.3. Combine to form a complete TM
Though the maximal target |bn(i)| on the output field can be obtained by each individual optimal incident field , we cannot get the maximal target |bn| when we set all optimal input patterns on the input field simultaneously. The reason might be that the maximal target |bn(i)| for each optimal input pattern has different phase, and they will have destructive interference on the output field. According to the theory that the TM is a linear matrix,[17] we find that the maximal target bn(i) will have a constant phase shift if we add the same constant phase shift to each element of optimal input pattern . As a result, we can search for the optimal phase shift of each optimal input pattern to achieve constructive interference on the target.
First, randomly selected phase shifts ( are added to each optimal input pattern. The number of randomly selected phase shift fields is m′ > 16log(16), and is all randomly selected phase shift fields. Then, we can get
where is the intensity of target for all input patterns. And we set as the phase of target for all input patterns. Then the optimal phase shift can be obtained by solving the following problem:
We also solve the above problem using the SDP method[13] and the PhaseCut algorithm,[16] such as Eq. (6)–(9). Then the optimal phase shift for each optimal input pattern is written as
Finally, the optimal input pattern for the whole input field can be written as
3. Experimental result
The experimental setup is shown in Fig. 3. We use a 532-nm semiconductor laser as the light source with 20-mW output power. The output beam is changed to vertical polarization state by the polarizer P1. Then the output beam from P1 is expended by the expending lens EL. The λ/2 wave plate and P2 are only used to change the intensity of light that illuminates on the DMD (Texas Instruments DLP 6500, 1920 × 1080 pixels). The wavefront of the reflected light is controlled by superpixel method. The lenses L1, L2, and the spatial filter SF make up a 4-f system. The +1 diffraction order of the reflected light can only pass the spatial filter SF. In the 4-f setup, the diameter of the output beam is half of that of the input beam. The modulated light illuminates on the scattering sample via a 40× objective Obj1. The scattering sample is the ground glass diffuser (Edmund #45-653). The output beam from the scattering sample is collected by a 10× objective Obj2 and illuminates on a CCD camera (AVT Manta G-031B, 656 × 492 pixels, pixel size 5.6 μm × 5.6 ×m).
We calculate a TM with 256 elements. First of all, we divide the TM into 4 × 4 parts and each part has 4 × 4 elements. In the same way, we divide the DMD into 4 × 4 parts and each part has 4 × 4 segments. Then each part of TM is measured in proper sequence as written in the Subsection 2.2. Finally, every part makes up a complete TM by a certain rule and we focused light to target though scattering media by computing optimal input pattern. The target intensity is recorded for each measurement in Fig. 4(a). The total number of measurements for each part is 45, so all parts require 720 measurements. The complete TM requires 45 measurements. Accordingly, the total number of measurements is 765. And the 766th intensity measurement could get the maximal target value for the whole optimal input pattern. The definition of the enhancement factor is η = Iopt/Iref,[7] where Iopt is the intensity of the target for the whole optimal input pattern, and Iref is the average intensity of the background. In Fig. 4(b), we show speckle near target before computing the optimal pattern. In Fig. 4(c), we show the maximal target intensity for the optimal pattern. The average of the enhancement factors obtained from ten experiments is 98.80.
Fig. 4. (color online) (a) The intensity at the target for 766 measurements on the CCD camera. (b) The target field before optimization. (c) The target field after optimization.
We also perform some experiments about different numbers of TM elements. We use the SDP method to calculate TM when N = 9, 16, 25, 36, and 49, and use SSDP method to calculate TM when N = 36, 81, 144, 256, and 400. The SSDP method and the SDP method are compared in four aspects. Due to the limitation of computer performance, we can only calculate the N = 49 TM with the SDP method. Each data point is the average of ten experiments. As shown in Fig. 5(a), we compared the two methods in enhancement factor. The enhancements obtained with the SDP method are 4.71 ± 1.63, 9.15 ± 1.91, 13.40 ± 2.07, 18.28 ± 1.75, and 25.02 ± 2.32, while the enhancements obtained with the SSDP method are 9.88 ± 2.5, 27.16 ± 5.6, 49.69 ± 1.8, 98.80 ± 6.9, and 139.89 ± 8.5. As can be seen, with the same segments number, the enhancement factor obtained with the SDP method is larger than that with the SSDP method. However, in the other three aspects, the performance of SDP method is not as good as that of SSDP. As shown in Fig. 5(b), the number of measurements of the SSDP method is less than that of the SDP method with the same number of segments. The number of measurements of the SDP method is 21, 46, 82, and 131, while the number of measurements of the SSDP method is 71, 199, 362, 766, and 1191. Figure 5(c) shows that the SSDP method costs much less calculation time than the SDP method. The calculation time of the SDP method is 2.3 s, 7.5 s, 55 s, 724 s, and 6359 s, while the calculation time of the SSDP method is 8 s, 18 s, 42 s, 174 s, and 1045 s. Figure 5(d) shows that the SSDP method consumes much less memory space than the SDP method. The SDP method consumes 4 Mb, 30 Mb, 321 Mb, 2153 Mb, and 29935 Mb of memory space, while the SSDP method consumes 4 Mb, 4 Mb, 30 Mb, 30 Mb, and 320 Mb. Although the enhancement factor obtained with the SDP method is larger than that of the SSDP method with the same number of segments, the SSDP method is more effective than the SDP method considering the limited computer resource.
Fig. 5. (color online) The SSDP method and the SDP method are compared in the aspects of (a) enhancement, (b) the number of measurements, (c) calculation time, and (d) memory space.
4. Conclusion
In this paper, we put forward a combined transmission matrix method to measure the monochromatic TM of scattering media for the first time, i.e., the SSDP method. The calculation speed of the SSDP method is investigated and compared with the SDP method in experiment. We find that the calculation speed of the SSDP method is much faster than that of the SDP by comparing the time consumed. In addition, the SSDP method requires less computer memory. The method also provides a new experimental idea to measure large TM of scattering medium. In our experiment, by combining the SSDP method and the superpixel method, we have realized rapid light focusing through the scattering medium. The phase of the input light field is modulated by a DMD (a pure amplitude modulation device) with the superpixel method. Using the measured TM, the optimal mask is computed. We believe that the SSDP method requires less calculation time when we use high-performance computer clusters. The SSDP method has shown the potential for imaging or light delivery through biological tissues.
